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ABSTRACT 

The head-on collision of equal sized drops is studied by full numerical simulations. 
The Navier-Stokes equations are solved for the fluid motion both inside and outside 
the drops using a front tracking/finite difference technique. The drops are accelerated 
toward each other by a body force that is turned off before the drops collide. When the 
drops collide, the fluid between them is pushed outward leaving a thin layer bounded 
by the drop surface. This layer gets progressively thinner as the drops continue to 
deform and in several of our calculations we artificially remove this double layer once 
it is thin enough, thus modeling rupture. If no rupture takes place, the drops always 
rebound, but if the film is ruptured the drops may coalesce permanently or coalesce 
temporarily and then split again. 

I. INTRODUCTION 

The dynamics of fluid drops is of considerable importance in a number of engineer- 
ing applications and natural processes. The combustion of fuel sprays, spray painting, 
various coating processes, as well as rain are only a few of the more common examples. 
While it is often possible to focus attention on the dynamic of a single drop and how 
it interacts with the surrounding flow, it is necessary to account for the interactions 


1 


between the drops and their collective effect on the flow when the number of drops per 
unit volume is high. The collision of two drops is an extreme case of two drop interac- 
tion and has been the topic of several investigations. The collision process generally 
involves large deformations and rupture of the interface separating the drops, and 
has not been amenable to detailed theoretical analysis. Previous studies are therefore 
mostly experimental, but sometimes supplemented by greatly simplified theoretical 
argument. Here, we present numerical simulations of the head-on collision of two 
drops, where the full Navier-Stokes equations are solved to give a detailed picture of 
the flow during collision. 

Previous investigations of droplet collision have been motivated by raindrop for- 
mation (Brazier-Smith, Jennings and Latham, 1972; Spengler and Gokhale, 1973, 
and others), by efforts to predict the phase distribution in agitated liquid-liquid dis- 
persions (Park and Blair, 1975), by concern about blade erosion due to dispersed 
liquid drops in low pressure turbines (Ryley and Bennett-Cowell, 1967) and by fuel 
spray behavior (Ashgriz and Givi, 1987). Recent experimental studies include those 
of Azhgriz and Poo (1990), and Jiang, Umemura and Law (1992) who show several 
sequences of photographs of the various mode of collision for both water and hydro- 
carbon drop. Drop collisions can generally be classified into four main categories: 
bouncing collision, where the drops collide and separate, retaining their identity; coa- 
lescence collision, where two drops become one; separation collision, where the drops 
temporarily become one but then break up again; and shattering collision, where the 
impact is so strong that the drops break up into several smaller drops. The form of 
the collision depends on the size of the drops, their relative velocities and the phys- 
ical properties of the fluids involved. In addition to head-on collisions, off-centered 
collisions (where the drops are displaced radially with respect to each other) is dis- 
cussed by both Azhgriz and Poo and Jiang et al. Other investigations may be found 
in Bradley and Stow (1978), Podvysotsky and Shraiber (1984) and Ashgriz and Givi 
(1987), for example. The major goal of these investigations has been to clarify the 
boundaries between the major collision categories and explain how they depend on 
the properties of the problem. As the large differences between water drops and hy- 
drocarbon drops, studied by Jiang et al., show, even the relatively simple case of head 
on collisions of equal size drops is still not fully understood. 

Previous theoretical investigations of drop behavior have almost all been concerned 
with the oscillations of a single drop. The linear oscillations of inviscid drops are well 
understood (see e.g. Lamb, 1932) and several authors have looked at nonlinear effects. 
Recent work includes analysis by Tsamopoulos and Brown (1983) and computations 
by Patzek, Benner, Basaran and Scriven (1991). The decay of linear oscillations due 
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to viscosity was analyzed in an approximate way by Lamb (1932) in the limit of small 
viscosity and a more detailed analysis was later carried out by Reid (1960), Miller and 
Scriven (1968) and others. Numerical investigations of viscous effects can be found in 
Foote (1973) who used the Marker And Cell (MAC) method to solve the full Navier 
Stokes equations, and Mansure and Lundgren (1988) who used a boundary integral 
method, modified to account for small viscous dissipation in an approximate way. 

The only simulation of drop collision that we are aware of is by Foote (1975) 
who followed the evolution of rebounding drops at low Reynolds and Weber number 
using the MAC method. The shapes computed by Foote have been compared with 
experimental observations by Bradley and Stow (1978) who found good agreement, 
but made the interesting observation that “this complicated treatment gives little 
insight into the physical processes involved.” Our simulations extend the work of 
Foote, increasing the Weber number and exploring what happens when the drops 
coalesce. 

The rest of the paper is laid out as follows: In section 2 we discuss briefly the 
numerical method which has been described in more detail elsewhere. Section 3 
contains our results and section 4 is devoted to discussions. In section 5 we summarize 
our results. Preliminary results have been presented at the 45 Annual Meeting of the 
Fluid Dynamics Division of the American Physical Society (Nobari and Tryggvason, 
1992). 

II. FORMULATION AND NUMERICAL METHOD 

The numerical technique used for the simulations presented in this paper is a front 
tracking method for multi-fluid flows developed by Unverdi (1990) and discussed by 
Unverdi and Tryggvason (1992 a,b). The actual code is an axisymmetric version of 
the method, described in Jan and Tryggvason (1993). Here we only briefly outline 
the procedure. 

The physical problem and the computational domain is sketched in Figure 1. The 
domain is axisymmetric and the drops are initially placed near each end of the domain. 
A force is applied to the drops to drive them together and turned off before the drops 
collide. Generally, the density and viscosity of the ambient fluid are much smaller 
than of the drop fluid and thus have only a small effect on the results. While it is 
therefore sufficient to solve only for the fluid motion inside the drop, here we solve 
for the motion everywhere, both inside and outside the drops. The Navier-Stokes 
equations are valid for both fluids, and a single set of equations can be written for 
the whole domain as long as the jump in viscosity and density is correctly accounted 
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for and surface tension is included. In component form these equations are 
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Here, v z and v r are the velocity components in the axial and radial direction, p is the 
pressure, and p and p are the discontinuous density and viscosity fields, respectively. 
F s is the surface tension force and f z is a body force. Notice that the surface tension 
force has been added as a delta function, only affecting the equations where the 
interface is. The detailed form of F s will be discussed below. The above equations 
are supplemented by the incompressibility conditions 
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which, when combined with the momentum equations leads to a non-separable elliptic 
equation for the pressure, and equations of state for the density and viscosity: 
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These last two equations simply state that density and viscosity within each fluid 
remains constant. 

Nondimensionalization gives a Weber and a Reynolds number defined by: 


We = 


p d dU 2 
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In addition, the density ratio p 0 l Pd and the viscosity ratio p 0 / Pd must be specified. 
Here, d is the drop diameter. The subscript d denotes the drop fluid and o the 
ambient fluid. When presenting our results we scale lengths by the initial diameter 
of the spherical drop and velocity by V = U/2, the speed of one drop before impact. 
To nondimensionalize time we have the choice of two inherent time scales: One is the 
advection time d/V of the drops before impact and the other is the natural oscillation 

time for the drop = (7r/4 )ypd^/a. While most of our results are presented using 
the advective time scale, in some cases the latter is the more natural one (as pointed 
out already by Foote, 1975). In a few cases we therefore replot our results using this 
time scale. 
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The force used to drive the drops together initially is taken as 

fz = A(p - p 0 )sign(z - z c ) 

so the force acts only on the drops. Here A is an adjustable constant and 2 C is midway 
between the drops. This force is turned off before the actual collision takes place. In 
most of our simulations the drops are initially put about one diameter apart (two 
diameters between their centers) and A is varied to give different collision velocities. 
To make comparison between various runs easier, we set time equal to zero when the 
centers of the drops are one diameter apart. If the drops were exactly spherical, they 
would touch at this instant. In our case, since the drops are moving in an other fluid, 
they have generally deformed slightly before impact and there is therefore a layer of 
ambient fluid between them at this time. 

To solve the Navier-Stokes equations we use a fixed, regular, staggered grid and 
discretize the momentum equations using a conservative, second order centered dif- 
ference scheme for the spatial variables and an explicit first order time integration 
method. We have used second order time integration in other problems and generally 
find little differences for relatively short simulations times as those of interest here. 
The effect does show up in long time simulations and is usually accompanied by a fail- 
ure to conserve mass. In the computations discussed here, mass is always conserved 
within a fraction of a percent. The interface is represented by separate computational 
points that are moved by interpolating their velocity from the grid. These points are 
connected to form a front that is used to keep the density and viscosity stratification 
sharp and to calculate surface tension forces. At each time step information must be 
passed between the front and the stationary grid. This is done by a method that has 
become known as the Immersed Boundary Method and is based on assigning the in- 
formation carried by the front to the nearest gridpoints. While this replaces the sharp 
interface by a slightly smoother grid interface, all numerical diffusion is eliminated 
since the grid-field is reconstructed at each step. 

The original Immersed Boundary Method was developed by Peskin and collabo- 
rators (see e.g. Peskin 1977) for homogeneous flows. The extension to multi-fluid 
flows includes a number of additional complications. The first is that density now 
depends on the position of the interface and has to be updated at each time step. 
There are several ways to do this but we use a variant of the method developed by 
Unverdi (1990) where the density jump at the interface is distributed onto the fixed 
grid to generate a grid-density-gradient field. The divergence of this field is equal to 
the Laplacian of the density field and the resulting Poisson equation can be solved 
efficiently by a Fast Poisson Solver. The particular attraction of this methods is that 
close interfaces can interact in a very natural way, since the grid-density-gradients 
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simply cancel. Therefore, when two interfaces come close together the full influence 
of the surface tension forces from both interfaces is included in the momentum equa- 
tions, but the mass of the fluids in the thin layer between the interfaces — which is very 
small — is neglected. A second complication is that the pressure equation now has a 
nonconstant coefficient (or is non-separable) since the density varies. This prevents 
the use of Fast Poisson Solvers based on Fourier Methods, or variants there of, and 
we have used a multigrid package, MUDPACK, from NCAR (see Adams, 1989, for a 
description) with slight modifications due to our staggered grid. 

The computation of the surface tension forces poses yet another difficulty. Gener- 
ally, curvature is very sensitive to minor irregularity in the interface shape and it is 
difficult to achieve accuracy and robustness at the same time. However, by computing 
the surface tension forces directly by 


F s = a [ ^ ds , 

JAs os 

where s is the tangent vector, we ensure that the net surface tension force is zero, or: 

j> crnnda = 0. 


Here, n is the outward normal and n is the mean curvature. This is important for 
long time simulations since even small errors can lead to a net force that moves the 
drop in an unphysical way. 

Lastly, contrary to previous computations with the Immersed Boundary Method, 
the interface deformes greatly in our simulations and it is necessary to add and delete 
computational elements during the course of the calculations. While this is a major 
task for fully three-dimensional simulations, here the interface is simply a line and 
such modifications are a simple matter. 

The method and the code has been tested in various ways, such as by extensive grid 
refinement studies, comparison with other published work and analytical solutions 
(for details see Jan, 1993 and Nobari, 1993). Generally, both analytical solutions and 
other simulations are limited to relatively simple cases. We include one test in figure 
2 where we compare the oscillations of a single drop with analytical predictions. Here 
a single drop is perturbed slightly by the fundamental mode. The drop oscillates 
and the amplitude of the fundamental mode is plotted in the figure. The oscillation 
period is close to what is predicted by Lamb (formula number 10 on page 475) with 
t CO mpute/ r <i = 1-03, and the decay compares also well with formula 12 on page 
641 in Lamb. The envelope for the oscillations, as computed by Lamb’s equation is 
plotted in figure 2. We have compared several cases and find, as expected, that as the 
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perturbation amplitude and the viscosity becomes smaller, fully resolved simulations 
give results in close agreement with the theoretical predictions. For large amplitude 
perturbations, the oscillation frequency is also well predicted by Lamb’s formula, if 
the diameter of a sphere of the same volume as the drop is used. 

m. RESULTS 

We have computed a large number of collisions and in this section we show several 
results. First we consider collisions where the interface between the drops is not 
ruptured, and then we discuss collisions with interface rupturing. 

A. Bouncing drops 

Figure 3 shows the collision of two drops, at several times. Here, We = 32, 
Re = 98, pd/ p 0 = 15, and pd/ Po = 350. Initially, a constant force acts on the drops to 
accelerate them toward each other. When the drops are about half a diameter apart, 
the force is turned off, but the drops have acquired enough momentum to continue 
toward each other and collide. As the drops come in contact, the fluid between them 
is squeezed away and the drops bulge out at the equator of the combined fluid mass. 
The bulk of the fluid continues to move forward and then outward to the rim of the 
drop — which is now more disk-like — thus resulting in an indentation in the middle. 
Surface tension eventually inhibits further outward motion of the rim and forces the 
fluid back toward the axis of symmetry. While kinetic energy is converted into surface 
tension energy during the initial deformation, the reversed motion converts surface 
tension energy back into kinetic energy and the drops rebound since the interface 
between the drops is not allowed to rupture. 

This calculation was done on a uniform grid with 64 by 256 meshes in the radial 
and axial direction, respectively. To show that this is an essentially fully converged 
solution we compare selected frames from the run in figure 3 with computations done 
on a coarser, 32 by 128, grid in figure 4. The most significant difference is that the 
coarsely resolved drops have moved slightly less apart than the well resolved ones, 
suggesting slightly larger loss of energy for low resolution. For all our simulations 
we have monitored the volume conservation of the drops (not explicitly enforced in 
the code) and found that even for collisions involving severe deformations the volume 
change is always less than a fraction of a percent. 

In figure 5 we show the time evolution of various quantities for the run in figure 
3. In (a) and (b) we show the position of the center of mass for the drop and the 
velocity of the center of mass, respectively. For reference, a horizontal line is drawn one 
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drop radius from the symmetry plane in (a), and if the drops remained spherical they 
would not touch once the center of mass is above this line. In reality the drops deform 
slightly before colliding and elongate after collision. The drops are set in motion by 
a constant force and the velocity in (b) therefore increases linearly. Once the force is 
turned off, the drops slow down slightly due to the drag from the surrounding fluid, 
and after collision the drops return to a steady motion with a velocity that is about 
half of their velocity before collision. 

In 5(c) we plot the nondimensional maximum diameter of the drop, and its de- 
formation, defined as the ratio of maximum diameter divided by the thickness of the 
drop on the symmetry axis. While the drops are squashed together during the colli- 
sion, they are elongated during rebounding. Notice that the radius does not start to 
increase until well into the collision and that the deformation curve is not symmetric 
about its peak, illustrating that the initial squashing is considerably slower than the 
subsequent recovery. The energy balance during the collision is shown in figure 5(d) 
where the kinetic energy of the drop and the surface tension energy, along with the 
total energy of the drop, is plotted versus time. Initially, only the kinetic energy 
increases as the drops are set in motion by the applied force field. When the force 
field is turned off, the energy decreases slightly due to viscous dissipation. During 
the collision the kinetic energy is converted into surface tension energy, which reaches 
maximum at the maximum drop deformation when the kinetic energy is nearly zero. 
Notice the rapid decay in the total energy during the initial stage of the collision when 
the drops become flatter; the “pause” in the dissipation at maximum deformation, 
and again the rapid dissipation when the drops recover their spherical shape. In the 
particular case shown here, the total energy dissipation is divided unequally between 
the initial deformation phase and the recovery phase, with larger dissipation taking 
place during the recovery stage. After the drops rebound, the surface and kinetic 
energy curves are slightly wavy, due to the oscillations of the drops and the total 
energy decreases at a rate comparable with the decay before collision. Comparing 5 
(a), (b), (c), and (d) we see that the maximum surface energy coincides with zero 
centroid velocity, but the maximum deformation occurs slightly later since the drops 
continue to become thinner in the center, even after the outward motion has stopped. 

Drops colliding with solid surfaces can cause extensive damage, and since the 
drops in figure 3 are of the same size, they behave as a drop colliding with a non- 
wetting, full-slip surface. Although real surfaces are not full-slip, it is likely that the 
differences are small at high Reynolds numbers. Figure 5e shows the force on the 
symmetry plane computed by integrating the pressure over the area where the drops 
touch. Notice that the maximum force coincides with the point where the centers of 
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mass are closest and that the maximum force is more than twice the average force. 
We have also estimated the average force by dividing the impulse needed to change 
the momentum of the drop by the collision time and find good agreement. 

To get additional insight into the collision process, we plot the velocity vectors 
inside the drop, as well as the streamlines for the whole flow field, at several times, in 
figure 6. The velocity vectors are to the right and (the mirror image of) the streamlines 
on the left. In the first frame the drops have collided and while the fluid in most of 
the drop is still moving forward with uniform velocity, the fluid in a small region near 
the collision plane, is moving outward. The forward motion of the drops has induced 
a circulation in the whole fluid domain leading to closed streamlines. In the outer 
fluid, near the drop surface there is a thin shear, visible as a “kink” in the streamlines. 
In the next frame, the region where the velocity is uniform and the streamlines are 
straight has nearly disappeared as more and more of the fluid is squeezed outward. 
Near the rim of the resulting disk the outward velocity eventually goes to zero, and 
in the third frame the outer rim is starting to flow inward, even though the middle 
of the disk is still getting thinner (the droplet never becomes completely stationary, 
thus the kinetic energy is never exactly zero). This reversed flow region continues to 
grow and in the fourth frame the flow is dominated by a large recirculation region of 
opposite circulation to the initial one. This development continues in the next two 
frames as the drops rebound. Since the flow near the walls of the domain is now 
toward the collision plane, a small reminder of the fluid with the original circulation 
accumulates near the outer walls. Notice that the flow field during recovery is not 
simply the reverse of the initial flow. While the drop was getting flatter, considerable 
amount of the drop fluid remained in uniform motion during a large fraction of the 
collision phase; during recovery the streamlines bend more uniformly. 

The pressure field inside the drops, at the same times as in figure 6, is plotted in 
figure 7. Because of finite resolution, the pressure is not exactly discontinuous across 
the interface, but changes smoothly over two to three grid spaces. For relatively fine 
resolution, as is the case here, this transition zone is thin. Initially, the pressure is 
nearly uniform within the drops, but as the drops collide and are brought to a halt, 
the pressure on the centerline, at the point of contact, increases. As the contact 
region increases the high pressure area moves to the rim of the disc, and at maximum 
deformation, when the drop is nearly stationary, the pressure is highest in the outer 
torus, where the curvature is highest. This high pressure drives the flow back during 
rebounding and as the drops separate the high pressure region is again on the contact 
plane. Here the drops are elongated during separation and the pressure is therefore 
highest near the ends where the curvature is highest. Notice that the vertical scale 
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in each frame is different. 

In laboratory experiments, bouncing is actually rather rare and drops generally 
coalesce when the film between them ruptures. The time of rupture depends on how 
fast fluid is drained out of the film and while we do not attempt to resolve the thin film, 
the computations give information that could be used to predict the drainage time. 
^From figure 3 and 4, it is clear that the area where the drops are in contact (separated 
only by a thin film) is well defined and not sensitive to the resolution. In figure 8a, we 
plot the radial position of the outer edge of this area versus time. In figure 8b we plot 
the pressure in the symmetry plane versus radius at several times. The earliest line, 
at t=0.2, corresponds to the second frame in figure 7. Here, the pressure is nearly 
uniform across the line of contact. As the drops become flatter the pressure on the 
centerline falls rapidly, but the pressure near the rim of the disk increases. This high 
pressure near the rim eventually forces the fluid back, again increasing the centerline 
pressure. As the drops rebound the centerline pressure becomes nearly uniform over 
the plane of contact and decreases. The radial velocity in the symmetry plane is 
plotted in figure 8c for the same times as in figure 8b. Initially, the velocity is nearly 
a linear function of radius, but then the fluid near the rim slows down and starts 
to flow back while the velocity near the center is still positive. Near the end of the 
collision, the velocities are again nearly linear, but now in the inward direction. This 
velocity is taken directly from the simulations in figure 3 and since the film is many 
times smaller than the grid spacing, this velocity is not affected by the properties 
of the film fluid. If the radial velocity in the film was different, one would expect a 
boundary layer inside the drop to match this velocity. 

Experimental observation suggest that the effect of the Reynolds number is small, 
once it is high enough. Although our Reynolds numbers are somewhat lower than 
those often encountered experimentally, we find a similar trend. In figure 9 we com- 
pare the results for a single Weber number and three Reynolds numbers. (The case 
shown in figure 3 has the same Weber number and is also included in the compar- 
isons made in figure 10 and 11.) Except for the very lowest Re, the solutions are 
quite similar. A more detailed comparison is given in figure 10 where the energies are 
plotted versus time. (Since the Weber numbers are not exactly identical, the curves 
do not coincide completely at time zero.) In all cases, the kinetic energy is reduced 
at nearly the same rate, but the rate of increase in surface tension energy is the same 
only for the two highest Reynolds numbers. The difference is due to dissipation as 
seen in the graph for the total energy. When the drops bounce back, the rate of 
recovery of kinetic energy is different for all Reynolds numbers, although the differ- 
ence is the smallest for the highest Reynolds numbers. Notice that for the highest 
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Reynolds numbers more energy is dissipated during the rebounding stage than during 
the initial impact. This is also seen in Figure 11a where the energy loss during initial 
impact (up to the maximum deformation) and the total energy loss are compared for 
the different Reynolds numbers. In figure lib the maximum radius and the maxi- 
mum deformation are plotted versus Reynolds number. While the maximum radius 
is relatively independent of Reynolds number once it is high enough, the deforma- 
tion continues to increase. Overall, the collision becomes relatively independent of 
Reynolds number for the highest values simulated and this is also reflected by the 
coefficient of restitution which is plotted in figure 11c. The restitution coefficient and 
the energy loss are computed when the distance between the drop center of mass is 
one diameter, since there is a small energy dissipation after the drops separate due 
to friction from the outer fluid. 

With collision at high Reynolds numbers becoming relatively independent of Reynolds 
number, the Weber number remains the main controlling parameter. Its influence on 
the collisions is examined in figure 12 where the drops are shown at several times for 
three different Weber numbers. In the top row the Weber number is smaller than 
in the computations in figure 3 and 4, but in the two lower rows the Weber num- 
bers are larger. There are obviously considerable differences. For the lowest Weber 
number the drops deform only slightly during the collision and return to nearly spher- 
ical shape immediately following separation. As the Weber number is increased the 
deformation increases considerably and the drops become greatly elongated as they 
separate. The time here is scaled by the initial velocity and on this timescale the 
collision takes longer for the higher Weber numbers. We have run the code at higher 
Weber numbers, but generally found it difficult to follow the computations throughout 
a complete bouncing due to instabilities in the thin film near the centerline. Whether 
this is a resolution problem or due to a physical instability has not been resolved. 
The question is most likely of marginal physical relevance since very thin films are 
likely to rupture for these high Weber numbers. 

The velocity of the center of mass is plotted in figure 13a. For the lowest Weber 
numbers the velocity changes smoothly from positive to negative, indicating a nearly 
constant deceleration of the center of mass. As the Weber number increases, the 
velocity decreases more rapidly, but the curve develops a kink at the point of maxi- 
mum deformation, where the velocity of the center of mass remains essentially zero 
as the drops become flatter. This “waiting” becomes longer as the Weber number 
increases and the final velocity of the drops after rebounding decreases due to the 
larger dissipation in the more deformed drops. The time in figure 13a is scaled using 
the velocity of the drops before collision. In figure 13b the time has been rescaled by 
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the oscillation frequency of the drop and in this units the total collision time is nearly 
constant! 

Figure 14 shows the force on the symmetry plane versus time. As the Weber num- 
ber increases, the drops become “softer” and the maximum lower. For the lowest 
Weber number the force has a single maximum, but for the higher Weber numbers 
there is a large maximum at the initial impact and another smaller one as the drops re- 
cover their shape and bounce back. The average force also decreases since the contact 
time increases and the net change of momentum during the collision becomes smaller 
since the final velocities are lower due to larger dissipation for larger deformation. 

In figure 15 we examine the energies as a function of time for the runs in figure 12 
using time units based on the oscillation frequency of a single drop (as in figure 13b). 
As the drops collide the kinetic energy (figure 15a) is reduced to nearly zero for all 
cases but the amount recovered depends strongly on the Weber number, with most 
energy dissipated for high Weber numbers where the deformations are large. This 
figure shows, as did figure 13b, that in the time units used here the time of collision 
is relatively constant for the higher Weber numbers. Furthermore, the post collision 
oscillations have nearly the same period — as expected. The surface tension energy, 
shown in figure 15b, actually has a smaller maximum for the higher Weber numbers, 
even though deformation is much larger, due to the increased dissipation. The total 
energy is shown in 15c. Here, it is obvious that the difference in energy losses is 
mostly due to different dissipation during the initial deformation. 

Figure 16 summarizes the results for different Weber numbers: As the Weber 
number increases, the drops deform more and the energy losses increase (16a), with 
nearly all the initial kinetic energy being dissipated at the highest We. The initial 
losses, up to maximum deformation, are about a third of the total losses for low We 
and increase to about half the losses for high We. As the deformation and energy 
dissipation increases, the restitution coefficient (16b) and the average collision force 
(16c) decrease. The collision time (16d), as measured in units based on the oscillation 
period of a single drop — and defined as the time from when the drops would first touch 
if they remained spherical until the time when the drops actually separate — decreases 
slightly at low Weber numbers and then remains relatively constant at higher We. 
This remarkable simple dependency of the collision time on We has been observed 
before (see e.g. Foote, 1975). For bouncing drops the collision time is, of course, of 
a critical importance, since it influences not only the total force exerted by the drop, 
but may also be important for mass and heat transfer. Furthermore, for coalescence 
to take place it is necessary that the collision takes sufficiently long time so that 
fluid can be drained from the film separating the drops. Translated into dimensional 
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variables, constant ^ co ili s i on / T i means, for example, that for a given fluid and drops 
size the collision time does not depend on the velocity of the drops. Low impact 
velocities (low We) will lead to small deformations, and large velocities (high We) to 
large deformations but the time in contact is the same. However, for the same fluids 
and same impact velocities, larger drops will have a longer contact time. Similarly, 
for the same size and impact velocities, drops with higher surface tension will bounce 
off each other faster than low surface tension drops. 

B. Coalescing drops 

In the above computations we have not ruptured the layer between the drops 
and therefore the drops can never coalesce. Real drops, however, generally coalesce 
(bouncing is actually somewhat rare) and the interface has to be ruptured for simula- 
tions of realistic collisions. Thin films usually rupture when their thickness becomes 
comparable with the intermolecular spacing (about 100-400 Angstrom, see, for ex- 
ample, Bradley and Stow 1978). We can not resolve the layer down to such a small 
scale, although the computations in figure 3-5 suggest that the large scale motion is 
well predicted and — in particular, and perhaps somewhat surprisingly — does not de- 
pend on the resolution of the layer. When the layer ruptures, however, the resulting 
change in the interface topology usually leads to dramatically different evolution from 
when the layer is not ruptured. The theory of film rupture between bubbles or drops 
is currently being developed (see e.g. Davis et al, 1989, and Yiantsios and Davis, 
1991), and while it appears possible that such a theory can be combined with full 
simulations, we take a more ad hoc approach here and rupture the interface at a pre- 
scribed time by removing surfaces which are very close. Such instantaneous change 
in topology is, of course, an approximation to what happens in reality. While the 
influence of molecular forces, where the actual rupture takes place, is confined to a 
small area, there generally follows an extremely rapid motion of the surrounding film 
where surface tension forces pull the remaining sheets and filament together, often 
leading to further rupture and the formation of small droplets. We ignore these rapid 
small scale processes entirely, also throwing away any small isolated drops that may 
be formed following the rupture. Modeling the rupture by a discontinuous change in 
the structure of the interface is therefore a little like modeling a shock wave by a dis- 
continuity. Although this “shock” is in time, rather than space, the analogy is made 
even more appropriate by the fact that usually the topology change is accompanied 
by a loss of surface and total energy. 

In figure 17 we show the same collision as in figure 3 and 4 where the interface 
is ruptured once the drops are close enough. In (a) the film is ruptured at time 
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0.4 by simply removing the double interface, leaving a single drop with indented 
waist. Surface tension pulls this indentation outward initially, but after the drop 
has reached its maximum deformation, surface tension pulls the waist inward and the 
drop elongates before starting to oscillate around the spherical equilibrium shape. The 
sensitivity of the evolution to the exact instance of rupture can be seen by comparing 
the frames in (a) to the frames in (b) where the interface is ruptured at a later time. 
The evolution is comparable to the previous case, but the maximum deformation 
is smaller. Figure 18 shows the evolution of the energies for the runs in figure 17, 
as well as the run in figure 3 when no rupturing takes place. As the interface is 
ruptured, considerable surface area disappears and there is therefore a discontinuous 
reduction in the surface energy (as well as the total energy). In reality this energy 
is dissipated when the ruptured film breaks into small drops or is stored as surface 
energy of these small drops, but here the film is simply removed. The kinetic energy 
is, of course, unchanged by the rupture, but its subsequent evolution is different than 
in the non-rupturing case. Notice that in (b) there is a larger energy loss and that 
the post-coalescence oscillations are smaller than in (a). 

We have repeated the computations in figure 12, where the Reynolds number is 
held constant (Re=100) and the Weber number varied, and ruptured the film between 
the drops at a predetermined nondimensional time (t=0.2). This early time was 
selected to minimize energy losses due to coalescence and since a well defined layer 
had formed at this stage so that removing it did not alter the total volume of the 
drop by any significant amount. For the We numbers simulated (up to 100) the drops 
coalesce permanently and figure 19a compares the maximum radius for these cases 
to the results where the drops bounce. When coalescence takes place, the maximum 
radius is larger. However, since some energy is lost when the thin film is removed, 
the maximum surface energy (19b) is smaller than for bouncing drops. 

Another simulation, for more energetic drops (Re=140, We=65) is shown in figure 
20 where we show the evolution following rupture for two different rupture times. In 
all cases the drops continue to become flatter, followed by a recovery that leads to a 
large elongation of the drop. For the first case where rupture is at an early time this 
elongation leads to a break up of the drop into two drops, but when the rupture is 
later this break up does not take place. In figure 21, the energies are plotted versus 
time. As the film is ruptured, there is a drop in surface energy and therefore total 
energy. Surface energy drops slightly following the rupture as the cusp left by the 
rupture is pulled back. The rate of decrease of kinetic energy is slowed, but not 
reversed, suggesting that considerable dissipation is taking place. As the combined 
drop continues to deform, surface energy increases again, reaching maximum where 
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the kinetic energy is minimum. Notice that the maximum is considerably later than 
when the interface is not ruptured. When the interface is ruptured earlier, the loss in 
energy is smaller, the maximum kinetic energy when the drop recovers its spherical 
shape is larger and subsequently, the surface energy at late time, when the kinetic 
energy has become nearly zero is slightly larger. This suggests that if the drops are 
allowed to coalesce earlier, a secondary separation will take place more easily. 

We have also conducted a few simulations at even higher Reynolds and Weber 
number. Figure 22 shows the evolution of the interface for Re=185 and We=115 
where the interface is ruptures at t=0.2. After coalescence and the initial formation 
of a flat “disk” the drops stretch apart, forming a chain of three nearly equal sized 
drops. Here, we have removed the filament connecting the drops after stretching, 
thus again modeled rupture. The size of the middle drop is considerably larger here 
than in figure 20. In experiments, several drops are often formed for more energetic 
collisions. 

IV. DISCUSSION 


In the modeling of droplet collisions the most basic question is what type of collision 
will result for a given set of external parameters. Most models proposed in the 
literature therefore try to predict the boundaries between the various collision modes. 
The simulations in the preceding section give detailed information about both the 
drop shape and the velocity field as a function of time and can help to validate the 
various hypotheses made in the construction of simple models. 

Both Ashgriz and Poo (1990) and Jiang et al (1992) present simple energy argu- 
ments to explain the outcome of drop collisions. The basic difference between these 
models is that Ashgriz and Poo neglect dissipative effects whereas Jiang et al include 
dissipation during deformation. For drops that coalesce, Jiang et al (1992) argue that 
the dissipation up to maximum deformation is independent of the viscosity of the 
fluid and that most of it takes place in a thin layer near the contact plane between 
the drops. 

^From figure 10 and 11 we see that while the collision becomes relatively indepen- 
dent of the Reynolds number as Re increases, the energy dissipation does not go to 
zero. Indeed, there seems to be some support for the assertion that the energy loss 
(particularly during the initial deformation) becomes independent of the Reynolds 
number. To examine this in a more detail, we plot the dissipation per unit volume 
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for selected times and three different Reynolds numbers in figure 23. The times were 
selected where the dissipation is high during the initial impact (t=0.2) and during 
rebound (t=1.2). The figure shows that the maximum dissipation does not take place 
in a thin layer near the stagnation point, as assumed by Jiang et al , but near the outer 
edge of the drop where the streamlines are turning outward. However, although the 
maximum dissipation is occurring in a different place than they assumed, the rest of 
their argument seems to be supported by the plot. While the contour plots for the 
highest Reynolds numbers, at t=0.2, are not identical, they are considerably closer to 
each other than to the plot for the lowest Reynolds number, thus suggesting some level 
of convergence. We note that this is actually a more stringent test than the argument 
of Jiang et al requires; here we are comparing the point-wise dissipation whereas their 
discussion is based on the integrated value. Similar trend is seen during the rebound 
stage (t=1.2) where the maximum dissipation takes place near the symmetry line 
away from the contact plane where the streamlines converge. Overall, the dissipation 
is not as localized as during the initial deformation and the differences between the 
plots for the highest Reynolds numbers are greater. Although energy dissipation 
during collision may become independent of Reynolds number for Re — > oo, we note 
that for coalescing drops, any excess energy must be dissipated by oscillations and 
the decay thus depend on Re. 

The dissipation of energy has a significant influence on the evolution of the drops 
after initial contact. In particular, large dissipation reduces the maximum deforma- 
tion. An upper bound on the maximum surface area can be easily determined (see 
e.g. Jiang et al , 1992): Since kinetic energy is converted into surface tension energy 
during collision, the surface area is maximum if no energy is lost and all the initial 
kinetic energy is converted into surface tension energy 

^ MdV 2 + aS 0 — CrSmax- 

Here, we ignore the outer fluid completely. Md is the mass of a single drop and S 0 
and Smax are the initial and maximum surface area, respectively. Assuming the drops 
to be spherical initially, and using the definition of the Weber number this can be 
written as 

_ trrSpV* _ We 
S 0 + 23cr47rr 2 + 48 ' 

This line is plotted in figure 24 for both the bouncing drops in figures 12-16 and 
the coalescing drops in figure 19. In both cases the maximum surface area is not 
achieved due to dissipation of energy. Since the interface is ruptured at a constant 
nondimensional time based on d/V (not oscillation period, tj , ) the drops are slightly 
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more deformed when the film is ruptured at higher Weber numbers and the difference 
between bouncing and coalescing drops therefore increases with We. In addition to 
our numerical results, we have also plotted data from Jiang et al (1992) in figure 24. 
The dotted line is a straight line fit to their data points. Overall there is a reasonable 
agreement (the data is, for example, bounded by our bouncing drops), but the slope 
of the experimental data is somewhat different than either of our curves. We expect 
that this is due to differences in the time of film rupture. At low Weber number, 
when the velocities are low, the time it takes to drain the film is likely to be long 
and losses due to rupture large. At higher We the opposite appears to hold. We note 
also that Jiang et al had to estimate the surface area from measurements of the drop 
radius, and some of the differences could be due to inaccuracies in this estimate. 

Computations at high Re and We require fine resolution and long computational 
time. We have therefore simulated only a few cases for reflective collisions, defined as 
when the drops separate following an initial coalescence. Using these few runs and 
experimental data from the literature we show, in figure 25, the boundaries between 
coalescence and reflective collisions in the Re- We plane. The crosses, that are con- 
nected by a solid line, are obtained from the data presented by Jiang et al (1992) and 
the line to the far right is from the high Reynolds number experiments of Ashgriz and 
Poo (1990). The circles represent our simulations. Open circles show a coalescence 
collision and filled circles stand for reflective collisions. In most cases the interface 
was ruptured at t=0.2. The experimental data does not extend to low Reynolds num- 
bers but our numerical data suggest — as one might expect — that reflective collisions 
do not take place at low Reynolds numbers. Although the comparison can only be 
qualitative — we do not, after all, have a physical model for the rupture time — the 
agreement is good where we have data and the numerical results suggest a natural 
extension of the experimental results to low Reynolds numbers. 

While the limited number of computations that we have done for reflective colli- 
sions does not allow us to draw general conclusions, the plot of the energies in figure 21 
suggest a relatively simple criteria for separation following initial coalescence: Com- 
paring the two graphs, we see that the surface tension energy during rebound exceeds 
that of two drops (the horizontal line) in (a) where the drops separate, but in (b) 
where the drops do not separate, the losses are sufficiently large so that surface ten- 
sion energy does not exceed that of two isolated drops. We therefore suspect that the 
drops will split if the losses due to coalescence and deformation are sufficiently small, 
or that 

2{]-M d V 2 + a So) - $ > 2aT 0 

Lj 

where S 0 is the surface area of a single spherical drop and is the total losses due 
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to both viscous dissipation and interface rupture. While the viscous losses are fully 
predicted by our computations, the losses due to rupture require accurate information 
about the time of rupture. 

V. CONCLUSIONS 

The computations of head-on collisions of two drops of equal size presented here 
are, in many ways, quite similar to those of Foote (1975) almost twenty years ago. 
Indeed, we have used many of the same diagnostics as those presented by Foote and 
been guided by his observations. We have, of course, been able to extend both the 
resolution of the calculations and the range of parameters examined. In particular, 
we have simulated much higher Weber numbers. We have also examined the effect 
of rupturing the interface between the drops at predetermined times and allowing 
coalescence. While these studies are not exhaustive and suffer from the lack of a 
physics based model for the rupture time, they show both that such a model is needed, 
and that the rest of the computations is under good control. While the details of the 
rupturing remain unresolved, the computations suggest that since the evolution is 
relatively insensitive to the resolution of the layer between the drops, the drainage 
process before rupture is primarily a one-way coupling in the sense that while the 
drop behavior affects the draining, the exact film behavior has minimal impact on 
the drop. The rupture time, on the other hand, is critical to the continuing evolution 
of the drop, and depends on how fast the film is drained. These observations suggest 
that a subgrid model, which takes in the pressure and velocity of the drop fluid and 
predicts the rupture time, which is the only information returned back to the drop 
simulations, would give a procedure that had a fully predictive capability. Such a 
subgrid model for the rupture, that is suitable for our approach, has been presented 
by Jacqmin and Foster (1993), but has not been incorporated into our code yet. We 
note that accurate prediction of the time of rupture requires careful tracking of the 
front and that numerical techniques that relay on grid based reconstruction of the 
interface (such as the Volume-of- Fluid method) are not able to predict the delay in 
rupture due to a finite drainage time. 

The simulations presented here are only a first step in a comprehensive numerical 
study of droplet collisions. Both a reliable rupture model and fully three dimensional 
simulations will be required before a complete insight and predictive capability are 
in place. Both extensions appear to be within sight. The subgrid model of Jacqmin 
and Foster (1993) was mentioned above; preliminary three-dimensional simulations 
are presented by Nobari and Tryggvason (1993). 
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Figure 1: The computational setup. The axisymmetric domain is bounded by full-slip 
walls and resolved by a regular grid. 
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Figure 2: Comparison with analytical predictions. Oscillations of a single drop. The 
computed oscillation frequency is 3.77. Linearized theory gives c j t h = 3.88 for completely 
inviscid drops. The rate of decay is also compared with the approximate theory of Lamb. 
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Figure 3: Collision of two drops. We = 32, Re = 98, pd/ p Q = 15, Pd/ Po = 350. The 
nondimensional time (scaled by the initial velocity and the drop diameter) is noted in 
each frame. The grid used here is 64 x 256 meshes. 



Figure 4: Resolution test. Selected frames from the computation in Figure 3 (left half) 
are compared with results obtained on twice as coarse grid (right half). The evolution 
on the coarser grid is slightly slower than on the finer grid. 
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Figure 5a: Diagnostics for the simulations in Figures 3 and 4. The position of the center 
of mass of the drops versus time. The horizontal line marks the distance when the centers 
are one diameters apart. 



Figure 5b: The velocity of the center of mass of the drops versus time. 
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Figure 5c: “Deformation” of the drop versus time (solid line) and maximum diameter 
(dashed line). 



Figure 5d: Kinetic energy, surface tension energy and total energy versus time 
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We = 32.60 



Figure 5e: The pressure force on the symmetry plane as a function of time. The horizontal 
line is the average force. 
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Figure 6: Velocity vectors inside the drops (left half) and streamlines (right half) for 
selected frames from the computations in figure 3. 
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Figure 7: The pressure for selected frames from the computations in figure 3. Notice 
that the vertical scale is different in each frame. The times are the same as in figure 6. 
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Figure 8a: The radial position of the outer edge of the contact plane between the drops 
as a function of time. 



Figure 8b: Pressure on the symmetry plane between the drops as a function of radius for 
several times. 
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Figure 8c: Radial velocity on the symmetry plane between the drops 
radius for several times. 
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Figure 9: Selected frames for the collision of two drops at We = 30, pd/ Po — 15, Pd/ Po — 
350, and different Reynolds numbers. The nondimensional time (based on initial velocity 
and drop diameter) is noted on each frame, (a) Re = 28. (b) Re = 58. (c) Re = 120. 
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Figure 10a: Kinetic energy for the runs in figure 9 and 3 versus time. 



Figure 10b: Surface energy for the runs in figure 9 and 3 versus time. 
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Figure 10c: Total energy for the runs in figure 9 and 3 versus time. 



Figure 11a: Diagnostics for the simulations in Figure 9 and 3. Loss of energy versus 
Reynolds number. The lower line shows the loss in total energy during first half of the 
collision (up to maximum deformation) and the top line shows the total loss during the 
collision. 
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Figure lib: Maximum deformation (solid line) and maximum radius (dashed line) versus 
Reynolds number. 



Figure 11c: Coefficient of restitution versus Reynolds number. 
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Figure 12: Selected frames for the collision of two drops at Re = 96, pd/ p 0 — 15, 
Pd/ p 0 = 350, and different We. The nondimensional time (based on initial velocity and 
drop diameter) is noted on each frame, (a) We = 13. (b) We = 66.9. (c) We = 112. 
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Figure 13a: Velocities versus time for the runs in Figures 12 and 3 versus nondimensional 
time. Time is nondimensionalized by the initial velocity of the drop and its diameter. 



Figure 13b: Time is nondimensionalized by the the fundamental period of the oscillations 
of a single drop. 
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Figure 14: Force on the symmetry plane versus time (scaled by the initial velocity of the 
drop and its diameter). 



Figure 15a: Kinetic energy for the runs in Figures 12 and 3 versus time. 
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Figure 15b: Surface energy for the runs in Figures 12 and 3 versus time. 



Figure 15c: Total energy for the runs in Figures 12 and 3 versus time. 
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Figure 16a: Diagnostics for the simulations in Figures 12 and 3. Loss of energy versus 
Weber number. The lower line shows the loss in total energy during first half of the 
collision (up to maximum deformation) and the top line shows the total loss during the 
collision. 



Figure 16b: Coefficient of restitution versus Weber number. 
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Figure 16 c: Average collision force versus Weber number. 



Figure 16 d: Time of collision in units of period of oscillation of a single drop 
Weber number. 


versus 
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(a) 



(b) 



Figure 17: The evolution following rupture of the interface separating the drops for the 
simulation in Figure 3. In both cases the drops coalesce permanently, (a) Rupture at 
t = 0.4. (b) Rupture at t = 0.6. 
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Figure 18a: The energy versus time for the simulations in figure 17a. 



Figure 18b: The energy versus time for the simulations in figure 17b. 
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Figure 19a: The evolution following rupture of the interface separating the drops at 
t = 0.2 for Re = 98, pdj p 0 = 15, and pd/ p 0 = 350 (same conditions as in figure 12). The 
maximum radius versus Weber number. 



Figure 19b: The maximum surface tension energy versus Weber number. 
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Figure 20: The evolution following rupture of the interface separating the drops for 
We = 65, Re = 140, Pd/Po = 15, and pd/po = 350. In (a) the drops eventually separate 
again, following initial coalescence, but in (b) the drops remain one. (a) Rupture at 
t = 0.2. (b) Rupture at t = 0.5. 
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Figure 21a: The energy versus time for the simulations in figure 19a as in Figure 20. 



Figure 21b: The energy versus time for the simulations in figure 19b as in Figure 20. 
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Figure 22: The evolution following rupture of the interface separating the drops for 
We = 115, Re = 185, p d / p 0 = 15, and p d /p 0 = 350. 
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Figure 23: Dissipation per unit volume for bouncing drops. t=0.2 for the left column 
and t=1.2 for the right column. Re=58 for the top row; Re=98 for the middle row; and 
Re=123 for the bottom row. 
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Figure 24: Maximum surface area. The top line is the theoretical prediction for no losses. 
The solid line is for bouncing drops and the dashed line is for drops that coalesce. The 
dash-dot line is a best fit to experimental data from Jiang et at. 



Figure 25: The boundaries between coalescing and separating collisions in the Re- We 
plane. Open circles are computations where the drops coalesced permanently, dark circles 
are computations where the drops separated again. The solid line is data from Jiang et 
at (1992) and the dashed line is an extrapolation based on the computational results. 
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